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The question of the nature of the dark matter in the Universe remains one of the 
most outstanding unsolved problems in basic science. One of the best motivated 
particle physics candidates is the lightest supersymmetric particle, assumed to 
be the lightest neutralino. We here describe DarkSUSY, an advanced numerical 
FORTRAN package for supersymmetric dark matter calculations which we release 
for public use. With the help of this package, the masses and compositions of 
various supersymmetric particles can be computed, for given input parameters of 
the minimal supersymmetric extension of the Standard Model (MSSM). For the 
lightest neutralino, the relic density is computed, using accurate methods which 
include the effects of resonances, pair production thresholds and coannihilations. 
Accelerator bounds are checked to identify viable dark matter candidates. Finally, 
detection rates are computed for a variety of detection methods, such as direct 
detection and indirect detection through antiprotons, gamma-rays and positrons 
from the Galactic halo or neutrinos from the center of the Earth or the Sun. 



1 Introduction 

One of the favourite candidates for the dark matter is a Weakly Interact- 
ing Massive Particle, a WIMP. In supersymmetric extensions of the standard 
model, the neutralino emerges as a natural WIMP candidate for the dark 
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matter of the universe. There is, however, not a unique way of extending 
the standard model with supersymmctry, but it is a general practice to use 
the simplest possible model, the minimal supersymmetric enlargement of the 
standard model (the MSSM), usually with some additional simplifying as- 
sumptions. 

Over several years, we have developed analytical and numerical FORTRAN 
tools for dealing with the sometimes quite complex calculations necessary to go 
from given input parameters in the MSSM to actual quantitative predictions 
of relic density in the universe of the neutralinos, and the direct and indirect 
detection rates. The program package, which we have named DarkSUSY, has 
now reached such a level of sophistication and maturity that we find the time 
appropriate for its public release. In the following sections, we briefly describe 
the different components of DarkSUSY, and refer the reader to our upcoming 
paper 1 for more details. 

For download of the latest version of DarkSUSY please visit the official 
DarkSUSY website, http : //www. physto . se/~edsjo/darksusy/. The version 
of DarkSUSY described here is 3.14.01-beta. 

2 Definition of the Supersymmetric model 

We work in the framework of the minimal supersymmetric extension of the 
standard model defined by, besides the particle content and gauge couplings 
required by supersymmetry, the superpotential and the soft supersymmetry- 
breaking potential. See 1,2 for details. 

It is obvious that the parameter space of the most general MSSM is huge, 
being specified by 124 a priori free parameters 3 . In DarkSUSY we reduce this 
set by assuming that the off-diagonal elements of the trilinear matrices A and 
the scalar mass matrices M are zero, and imposing CP conservation (except 
in the CKM matrix). 

We then calculate all masses and most of the vertices entering the Feyn- 
man rules, which are all available to the user. We include several options for 
the loop corrections to the Higgs masses. 4 ' 5,6 ' 7 We also include loop corrections 
to the neutralino and chargino masses. 8 ' 9 

3 Accelerator bounds 

Accelerator bounds can be checked by a call to a subroutine. By modifying 
an option, the user can impose bounds as of different moments in time. The 
default option in version 3.14.01-beta adopts the 2000 limits by the Particle 
Data Group 10 modified slightly for the Higgs masses. The user is also free 
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to use his own routine to check for experimental bounds, in which case he or 
she would only need to provide an interface to DarkSUSY. 



4 Calculation of the relic density 

In DarkSUSY, we calculate the relic density by using the full cross section, 
including all resonances and thresholds and solve the Boltzmann equation 
numerically with the method given in n < 12 . 

When other supersymmetric particles are close in mass to the lightest 
neutralino they will also be present at the time of freeze-out in the early Uni- 
verse. When this happens so called coannihilations can take place between 
all these supersymmetric particles. We include the coannihilation processes 
between all charginos and neutralinos lighter than f co m x . The mass frac- 
tion parameter f co is by default set to 2.1 or 1.4 depending on how the relic 
density routines are called (high accuracy or fast calculation), but can be set 
to anything by the user. We have not included coannihilations with squarks 
and staus which occurs more accidentally than the in many cases unavoidable 
mass degeneracy between the lightest neutralinos and the lightest chargino. 

We have included all two-body final states that occur at tree level, both 
for neutralino-ncutralino, neutralino-chargino and chargino-chargino annihi- 
lations. Annihilation to gg, 77 and Z7 that occur at the 1-loop level are also 
included. 



5 Detection rates 

The different detection rates for neutralino dark matter have been calculated 
by many authors in the past. We will here only give a brief review about 
what is included in DarkSUSY, and which calculations they are based on. For 
a more extensive list of references, we refer to 1 . 



5.1 Halo models 

Currently implemented in DarkSUSY is the spherical family of halo profiles 
p(r) oc l/((r/a) 7 [1 + (r/a) a ]^ i)/ a ^j where e.g. the Navarro, Frenk and 
White profile 13 is given by (a, (3, 7) = (1,3,1) and the isothermal sphere 
is given by (a, (3, 7) = (2,2,0). The velocity distribution is assumed to be a 
standard isotropic gaussian distribution. 
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5.2 Direct detection 



These routines calculate the spin-dependent and spin-independent scattering 
cross sections on protons and neutrons assuming the quark contributions to 
the nucleon spin from 14 . The older set of data from 15 is also available as an 
option, but the user can set their own values if they wish. 

5.3 Monte Carlo simulations 

In several of the indirect detection processes below we need to evaluate the 
yield of different particles per neutralino annihilation. The hadronization 
and/or decay of the annihilation products are simulated with Pythia 16 and 
the results are tabulated. These tables are then used by DarkSUSY. 

5.4 Neutrinos from the Sun and Earth 

Neutralinos can accumulate in the Earth and the Sun where they can anni- 
hilate pair-wise producing high energy muon neutrinos. The branching ratios 
for different annihilation channels are calculated and the Pythia simulations 
are used to evaluate the yield of neutrinos. Neutrino interactions in the Sun 
as well as the charged current neutrino-nucleon interaction near the detector 
are also simulated with Pythia. 

There are routines to calculate a) the neutrino flux, b) the neutrino-to- 
muon conversion rate and c) the neutrino- induced muon flux either differential 
in energy and angle or integrated within an angular cone and above a given 
threshold. The new population of neutralinos in the solar system (arising 
from neutralinos that have scattered in the outskirts of the Sun) as described 
in 17 > 18 > 19 can optionally be included as well. 

5.5 Antiprotons from halo annihilations 

Neutralinos can also annihilate in the Milky Way halo producing e.g. antipro- 
tons. These propagate in the galaxy before reaching us. We have implemented 
the propagation method described in 20 . Optionally, the antiproton fluxes can 
also be solar modulated with the spherically symmetric model of 21 . There are 
also other propagation models 22 available as options. The antiproton fluxes 
are given differential in energy. 

5. 6 Positrons from halo annihilations 

In neutralino annihilations in the halo we can also produce positrons. The 
flux of positrons is calculated with the propagation model in 23 (with two 
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choices of the energy dependence of the diffusion constant). The model in 24 
can also be used as an option. The positron fluxes are given differential in 
energy. 

5. 7 Gamma rays from halo annihilations 

We can also produce gamma rays from the halo annihilations. These arc cither 
monochromatic, produced from f-loop annihilation 25 into 77 and 26 Z7, or 
with a continuous energy spectrum, produced from ir° decays in quark jets 27 . 

The flux of gamma rays can be obtained in any given direction on the sky 
for the user's choice of (a, f3, 7) in the halo profile. There arc also routines 
to average the flux over a chosen angular resolution. The continuous gamma 
rays use Pythia simulations to calculate the gamma ray flux (differential in 
energy, or integrated above an energy threshold). 

5. 8 Neutrinos from halo annihilations 

We can also produce neutrinos from neutralino annihilations in the halo. Al- 
though the fluxes are small, there are routines to calculate a) the neutrino flux, 
b) the neutrino-to-muon conversion rate and c) the neutrino-induced muon 
flux either differential in energy and angle or integrated within an angular 
cone and above a given energy threshold. 

6 Conclusions 

Over the years we have developed this numerical package, DarkSUSY, for 
neutralino dark matter calculations in the Minimal Supersymmetric Stan- 
dard Model, MSSM. The package is now publically released and available 
for download from http://www.physto.se/~edsjo/darksusy. We have here 
presented an overview of what the program can do, and refer the reader to the 
upcoming paper , where the details will be given. A test program, provided 
with the distribution, shows in more detail how DarkSUSY is used. 
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